Analysis for SSFA project: Three factor model

Table of contents

1. Used packages
2. Global settings
3. Load data
4. Explore data
5. Model specification
6. Inference
   1. epLsar
   2. R²
   3. Asfc
   4. Smfc
   5. HAsfc9
   6. HAsfc81
7. Summary

Used packages

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import pickle
import arviz as az
import pymc3 as pm

WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

from matplotlib.colors import to_rgb

import scipy.stats as stats 

from IPython.display import display

import matplotlib as mpl

%load_ext autoreload
%autoreload 2

import plotting_lib

Global settings

Output

writeOut = True
outPathPlots = "../plots/statistical_model_three_factors/"
outPathData = "../derived_data/statistical_model_three_factors/"

Plotting

widthMM = 190 
widthInch = widthMM / 25.4
ratio = 0.66666
heigthInch = ratio*widthInch

SMALL_SIZE = 8
MEDIUM_SIZE = 10
BIGGER_SIZE = 12

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title
sns.set_style("ticks")

dpi = 300

sizes = [SMALL_SIZE,MEDIUM_SIZE,BIGGER_SIZE]

Computing

numSamples = 1000
numCores = 10
numTune = 1000
numPredSamples = 2000
random_seed=36514535
target_accept = 0.99

Load data

datafile = "../derived_data/preprocessing/preprocessed.dat"

with open(datafile, "rb") as f:
    x1,x2,x3,df,dataZ,dictMeanStd,dictTreatment,dictSoftware = pickle.load(f)    

Show that everything is correct:

display(pd.DataFrame.from_dict({'x1':x1,'x2':x2,'x3':x3}))

	x1	x2	x3
0	0	5	116
1	1	5	116
2	0	5	117
3	1	5	117
4	0	5	118
...	...	...	...
275	1	9	51
276	0	9	52
277	1	9	52
278	0	9	53
279	1	9	53

280 rows × 3 columns

x1 indicates the software used, x2 indicates the treatment applied and x3 the sample.

for surfaceParam,(mean,std) in dictMeanStd.items():
    print("Surface parameter {} has mean {} and standard deviation {}".format(surfaceParam,mean,std))

Surface parameter epLsar has mean 0.0028221512571428575 and standard deviation 0.0019000504797019124
Surface parameter R² has mean 0.998765273042857 and standard deviation 0.0015328558023807836
Surface parameter Asfc has mean 18.13129049385357 and standard deviation 16.348381888991312
Surface parameter Smfc has mean 4.938492967075 and standard deviation 38.106353908569346
Surface parameter HAsfc9 has mean 0.3085220006979256 and standard deviation 0.2211140516395699
Surface parameter HAsfc81 has mean 0.5639598041398011 and standard deviation 0.40668126467296095

for k,v in dictTreatment.items():
    print("Number {} encodes treatment {}".format(k,v))

Number 5 encodes treatment Dry Bamboo
Number 6 encodes treatment Dry Grass
Number 7 encodes treatment Dry Lucerne
Number 0 encodes treatment BrushDirt
Number 1 encodes treatment BrushNoDirt
Number 4 encodes treatment Control
Number 10 encodes treatment RubDirt
Number 2 encodes treatment Clover
Number 3 encodes treatment Clover+Dust
Number 8 encodes treatment Grass
Number 9 encodes treatment Grass+Dust

for k,v in dictSoftware.items():
    print("Number {} encodes software {}".format(k,v))

Number 0 encodes software ConfoMap
Number 1 encodes software Toothfrax

display(dataZ)

	index	TreatmentNumber	SoftwareNumber	DatasetNumber	NameNumber	epLsar_z	R²_z	Asfc_z	Smfc_z	HAsfc9_z	HAsfc81_z
0	0	5	0	0	116	0.839414	-0.017104	-0.448128	-0.120476	-0.806963	-0.512247
1	1	5	1	0	116	0.999368	0.518462	-0.477757	-0.126469	-0.782634	-0.497016
2	2	5	0	0	117	1.601888	-0.509024	-0.276513	-0.119325	-0.584158	-0.662886
3	3	5	1	0	117	1.596720	0.457791	-0.301685	-0.126469	-0.629947	-0.744422
4	4	5	0	0	118	1.168099	-0.221668	-0.393502	-0.121498	-0.269712	-0.370958
...	...	...	...	...	...	...	...	...	...	...	...
275	275	9	1	2	51	0.843056	0.387986	-0.997743	-0.092631	2.388080	1.346868
276	276	9	0	2	52	0.305544	-0.791837	-0.967607	-0.104937	2.014963	2.573677
277	277	9	1	2	52	0.166758	0.635237	-0.940337	-0.126098	2.926894	3.117333
278	278	9	0	2	53	-0.843412	0.042974	-1.022523	-0.085082	0.280534	0.543577
279	279	9	1	2	53	-1.115313	0.712218	-1.021455	-0.118879	0.169999	0.678630

280 rows × 11 columns

display(df)

	Dataset	Name	Software	Diet	Treatment	Before.after	epLsar	R²	Asfc	Smfc	HAsfc9	HAsfc81	NewEplsar	TreatmentNumber	SoftwareNumber	DatasetNumber	NameNumber
0	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_1	ConfoMap	Dry Bamboo	Dry Bamboo	NaN	0.004417	0.998739	10.805118	0.347586	0.130091	0.355639	0.019460	5	0	0	116
1	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_1	Toothfrax	Dry Bamboo	Dry Bamboo	NaN	0.004721	0.999560	10.320730	0.119219	0.135471	0.361833	NaN	5	1	0	116
2	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_2	ConfoMap	Dry Bamboo	Dry Bamboo	NaN	0.005866	0.997985	13.610750	0.391436	0.179356	0.294377	0.020079	5	0	0	117
3	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_2	Toothfrax	Dry Bamboo	Dry Bamboo	NaN	0.005856	0.999467	13.199232	0.119219	0.169232	0.261217	NaN	5	1	0	117
4	GuineaPigs	capor_2CC6B1_txP4_#1_1_100xL_3	ConfoMap	Dry Bamboo	Dry Bamboo	NaN	0.005042	0.998425	11.698166	0.308648	0.248885	0.413098	0.019722	5	0	0	118
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
275	Sheeps	L8-Ovis-90730-lm2sin-a	Toothfrax	Grass+Dust	Grass+Dust	NaN	0.004424	0.999360	1.819802	1.408678	0.836560	1.111706	NaN	9	1	2	51
276	Sheeps	L8-Ovis-90764-lm2sin-a	ConfoMap	Grass+Dust	Grass+Dust	NaN	0.003403	0.997552	2.312486	0.939718	0.754059	1.610626	0.018978	9	0	2	52
277	Sheeps	L8-Ovis-90764-lm2sin-a	Toothfrax	Grass+Dust	Grass+Dust	NaN	0.003139	0.999739	2.758297	0.133366	0.955699	1.831721	NaN	9	1	2	52
278	Sheeps	L8-Ovis-90814-lm2sin-a	ConfoMap	Grass+Dust	Grass+Dust	NaN	0.001220	0.998831	1.414701	1.696316	0.370552	0.785022	0.017498	9	0	2	53
279	Sheeps	L8-Ovis-90814-lm2sin-a	Toothfrax	Grass+Dust	Grass+Dust	NaN	0.000703	0.999857	1.432148	0.408433	0.346111	0.839946	NaN	9	1	2	53

280 rows × 17 columns

Contrasts

Prepare dicts for contrasts that are of interest later.

x1contrast_dict = {'{}_{}'.format(dictSoftware[0],dictSoftware[1]):[1,-1]}
x1contrast_dict

{'ConfoMap_Toothfrax': [1, -1]}

Explore data

Compute raw differences

dfRawDiff = dataZ.groupby(['DatasetNumber','TreatmentNumber','NameNumber']).agg(list).applymap(lambda l : l[0]-l[1]).reset_index()

dfRawDiff

	DatasetNumber	TreatmentNumber	NameNumber	index	SoftwareNumber	epLsar_z	R²_z	Asfc_z	Smfc_z	HAsfc9_z	HAsfc81_z
0	0	5	116	-1	-1	-0.159954	-0.535566	0.029629	0.005993	-0.024329	-0.015230
1	0	5	117	-1	-1	0.005168	-0.966815	0.025172	0.007144	0.045788	0.081537
2	0	5	118	-1	-1	0.139785	-0.550288	0.039022	0.004971	0.142403	0.139025
3	0	5	119	-1	-1	0.370716	-0.763882	0.041183	0.007772	0.113397	-0.223817
4	0	5	120	-1	-1	0.138629	-0.150835	0.074677	0.007772	0.179692	0.075084
...	...	...	...	...	...	...	...	...	...	...	...
135	2	9	49	-1	-1	0.576471	-0.505284	-0.008495	0.004560	0.011585	0.007728
136	2	9	50	-1	-1	1.728006	-1.316577	-0.005035	0.005599	0.579285	0.167894
137	2	9	51	-1	-1	-2.151764	-2.659380	-0.008650	0.510897	0.095887	0.165743
138	2	9	52	-1	-1	0.138786	-1.427074	-0.027269	0.021161	-0.911931	-0.543656
139	2	9	53	-1	-1	0.271901	-0.669244	-0.001067	0.033797	0.110535	-0.135053

140 rows × 11 columns

Show difference between software in raw data

variablesList = ['epLsar','R²','Asfc','Smfc','HAsfc9','HAsfc81']

for var in variablesList:
    fig, axes = plt.subplots(1, 2,figsize=(12,6))
    fig.suptitle('{}'.format(var))

    
    sns.stripplot(data=dataZ,x='NameNumber',y='{}_z'.format(var),hue='SoftwareNumber',ax=axes[0]);
    axes[0].set_title('Raw data')

    sns.histplot(data=dfRawDiff,x='{}_z'.format(var),ax=axes[1]);
    axes[1].set_title('Aggegrated difference by software'.format(var))

    plt.show()

Model specification

class ThreeFactorModel(pm.Model):
    
    """
    Compute params of priors and hyperpriors.
    """
    def getParams(self,x1,x2,x3,y):
        # get lengths
        Nx1Lvl = np.unique(x1).size
        Nx2Lvl = np.unique(x2).size
        Nx3Lvl = np.unique(x3).size        
        
        dims = (Nx1Lvl, Nx2Lvl, Nx3Lvl)
        
        ### get standard deviations
        
        # convert to pandas dataframe to use their logic
        df = pd.DataFrame.from_dict({'x1':x1,'x2':x2,'x3':x3,'y':y})
        
        s1 = df.groupby('x1').std()['y'].max()
        s2 = df.groupby('x2').std()['y'].max()
        s3 = df.groupby('x3').std()['y'].max()
        
        stdSingle = (s1, s2, s3)
        
        return (dims, stdSingle)
    
    def printParams(self,x1,x2,x3,y):
        dims, stdSingle = self.getParams(x1,x2,x3,y)
        Nx1Lvl, Nx2Lvl, Nx3Lvl = dims
        s1, s2, s3 = stdSingle        
        
        print("The number of levels of the x variables are {}".format(dims))
        print("The standard deviations used for the beta priors are {}".format(stdSingle))
            
    def __init__(self,name,x1,x2,x3,y,model=None):
        
        # call super's init first, passing model and name
        super().__init__(name, model)
        
        # get parameter of hyperpriors
        dims, stdSingle = self.getParams(x1,x2,x3,y)
        Nx1Lvl, Nx2Lvl, Nx3Lvl = dims
        s1, s2, s3 = stdSingle        
        
        ### hyperpriors ### 
        # observation hyperpriors
        lamY = 1/30.
        muGamma = 0.5
        sigmaGamma = 2.
        
        # prediction hyperpriors
        sigma0 = pm.HalfNormal('sigma0',sd=1)
        sigma1 = pm.HalfNormal('sigma1',sd=s1, shape=Nx1Lvl)
        sigma2 = pm.HalfNormal('sigma2',sd=s2, shape=Nx2Lvl)
        sigma3 = pm.HalfNormal('sigma3',sd=s3, shape=Nx3Lvl)
                
        mu_b0 = pm.Normal('mu_b0', mu=0., sd=1)
        mu_b1 = pm.Normal('mu_b1', mu=0., sd=1, shape=Nx1Lvl)
        mu_b2 = pm.Normal('mu_b2', mu=0., sd=1, shape=Nx2Lvl)
        mu_b3 = pm.Normal('mu_b3', mu=0., sd=1, shape=Nx3Lvl)       
                                       
        ### priors ### 
        # observation priors        
        nuY = pm.Exponential('nuY',lam=lamY)
        sigmaY = pm.Gamma('sigmaY',mu=muGamma, sigma=sigmaGamma)
        
        # prediction priors
        b0_dist = pm.Normal('b0_dist', mu=0, sd=1)
        b0 = pm.Deterministic("b0", mu_b0 + b0_dist * sigma0)
       
        b1_dist = pm.Normal('b1_dist', mu=0, sd=1)
        b1 = pm.Deterministic("b1", mu_b1 + b1_dist * sigma1)
        
        b2_dist = pm.Normal('b2_dist', mu=0, sd=1)
        b2 = pm.Deterministic("b2", mu_b2 + b2_dist * sigma2)
        
        b3_dist = pm.Normal('b3_dist', mu=0, sd=1)
        b3 = pm.Deterministic("b3", mu_b3 + b3_dist * sigma3)
             
        #### prediction ###      
        mu = pm.Deterministic('mu',b0 + b1[x1] + b2[x2] +  b3[x3])

                                
        ### observation ### 
        y = pm.StudentT('y',nu = nuY, mu=mu, sd=sigmaY, observed=y)

Inference

epLsar

with pm.Model() as model:
    epLsarModel = ThreeFactorModel('epLsar',x1,x2,x3,dataZ.epLsar_z.values)

Verify model settings

epLsarModel.printParams(x1,x2,x3,dataZ.epLsar_z.values)

The number of levels of the x variables are (2, 11, 140)
The standard deviations used for the beta priors are (1.0034852973454325, 1.2901096312620055, 1.5215268338889083)

try:
    graph_epLsar = pm.model_to_graphviz(epLsarModel)    
except:
    graph_epLsar = "Could not make graph"
graph_epLsar

Check prior choice

with epLsarModel as model:
    prior_pred_epLsar = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')

Prior choice is as intended: Broad over the data range.

Sampling

with epLsarModel as model:
    trace_epLsar = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
    #fit_epLsar = pm.fit(random_seed=random_seed)
    #trace_epLsar = fit_epLsar.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [epLsar_b3_dist, epLsar_b2_dist, epLsar_b1_dist, epLsar_b0_dist, epLsar_sigmaY, epLsar_nuY, epLsar_mu_b3, epLsar_mu_b2, epLsar_mu_b1, epLsar_mu_b0, epLsar_sigma3, epLsar_sigma2, epLsar_sigma1, epLsar_sigma0]

100.00% [20000/20000 18:55<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 1_000 draw iterations (10_000 + 10_000 draws total) took 1137 seconds.
The number of effective samples is smaller than 10% for some parameters.

with epLsarModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('epLsar'), 'wb') as buff:
            pickle.dump({'model':epLsarModel, 'trace': trace_epLsar}, buff)

Check sampling

with epLsarModel as model:
    dataTrace_epLsar = az.from_pymc3(trace=trace_epLsar)

pm.summary(dataTrace_epLsar,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
epLsar_mu_b0	-0.06	0.82	-1.67	1.48	0.01	0.01	8287.0	5027.0	8296.0	7497.0	1.0
epLsar_mu_b1[0]	0.00	0.70	-1.34	1.39	0.01	0.01	8365.0	5752.0	8342.0	7640.0	1.0
epLsar_mu_b1[1]	-0.06	0.69	-1.37	1.32	0.01	0.01	8990.0	5503.0	8990.0	7313.0	1.0
epLsar_mu_b2[0]	-0.26	0.53	-1.32	0.75	0.01	0.00	4827.0	4811.0	4780.0	6328.0	1.0
epLsar_mu_b2[1]	-0.50	0.53	-1.49	0.57	0.01	0.01	4168.0	4149.0	4156.0	5540.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
epLsar_mu[275]	0.43	0.76	-1.42	1.13	0.03	0.02	710.0	710.0	1342.0	2187.0	1.0
epLsar_mu[276]	0.27	0.12	0.04	0.50	0.00	0.00	11139.0	8611.0	10983.0	7449.0	1.0
epLsar_mu[277]	0.21	0.12	-0.03	0.43	0.00	0.00	11102.0	7679.0	10946.0	7539.0	1.0
epLsar_mu[278]	-0.92	0.13	-1.17	-0.65	0.00	0.00	10940.0	10940.0	11017.0	7556.0	1.0
epLsar_mu[279]	-0.98	0.13	-1.25	-0.73	0.00	0.00	10936.0	10936.0	11000.0	7728.0	1.0

748 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,dataTrace_epLsar,'epLsar')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with epLsarModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_epLsar,'epLsar')

with epLsarModel as model:
    plotting_lib.pm.energyplot(trace_epLsar)

Posterior predictive distribution

with epLsarModel as model:
    posterior_pred_epLsar = pm.sample_posterior_predictive(trace_epLsar,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_epLsar,posterior_pred_epLsar,dataZ.epLsar_z.values,'epLsar')

Compare prior and posterior for model parameters

with epLsarModel as model:
    pm_data_epLsar = az.from_pymc3(trace=trace_epLsar,prior=prior_pred_epLsar,posterior_predictive=posterior_pred_epLsar)

arviz.data.io_pymc3 - WARNING - posterior predictive variable epLsar_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_epLsar,'epLsar')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_epLsar,'epLsar')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_epLsar,'epLsar')

R²

with pm.Model() as model:
    RsquaredModel = ThreeFactorModel('R²',x1,x2,x3,dataZ["R²_z"].values)

Verify model settings

RsquaredModel.printParams(x1,x2,x3,dataZ["R²_z"].values)

The number of levels of the x variables are (2, 11, 140)
The standard deviations used for the beta priors are (1.0112351887037512, 2.890846695679254, 6.372712717079752)

pm.model_to_graphviz(RsquaredModel)

Check prior choice

with RsquaredModel as model:
    prior_pred_Rsquared = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,dataZ["R²_z"].values,'R²')

Sampling

with RsquaredModel as model:
    trace_Rsquared = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
    #fit_Rsquared = pm.fit(random_seed=random_seed)
    #trace_Rsquared = fit_Rsquared.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [R²_b3_dist, R²_b2_dist, R²_b1_dist, R²_b0_dist, R²_sigmaY, R²_nuY, R²_mu_b3, R²_mu_b2, R²_mu_b1, R²_mu_b0, R²_sigma3, R²_sigma2, R²_sigma1, R²_sigma0]

100.00% [20000/20000 10:34<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 1_000 draw iterations (10_000 + 10_000 draws total) took 636 seconds.
The number of effective samples is smaller than 10% for some parameters.

with RsquaredModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('R²'), 'wb') as buff:
            pickle.dump({'model': RsquaredModel, 'trace': trace_Rsquared}, buff)

Check sampling

with RsquaredModel as model:
    dataTrace_Rsquared = az.from_pymc3(trace=trace_Rsquared)

pm.summary(dataTrace_Rsquared,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
R²_mu_b0	0.02	0.80	-1.54	1.60	0.01	0.01	7488.0	6433.0	7488.0	7567.0	1.00
R²_mu_b1[0]	-0.21	0.69	-1.60	1.12	0.01	0.01	6630.0	6234.0	6640.0	7314.0	1.00
R²_mu_b1[1]	0.25	0.69	-1.09	1.59	0.01	0.01	6754.0	6464.0	6756.0	6869.0	1.00
R²_mu_b2[0]	0.02	0.52	-0.97	1.06	0.01	0.01	2703.0	2703.0	2701.0	5022.0	1.00
R²_mu_b2[1]	0.02	0.52	-1.02	1.03	0.01	0.01	2798.0	2798.0	2787.0	4860.0	1.00
...	...	...	...	...	...	...	...	...	...	...	...
R²_mu[275]	0.04	0.67	-1.74	0.77	0.02	0.01	1631.0	1535.0	3108.0	3027.0	1.01
R²_mu[276]	-0.32	0.38	-1.00	0.32	0.00	0.00	8948.0	6818.0	10701.0	7802.0	1.00
R²_mu[277]	0.22	0.38	-0.46	0.86	0.00	0.00	8920.0	6761.0	10664.0	8104.0	1.00
R²_mu[278]	0.10	0.15	-0.20	0.42	0.00	0.00	10768.0	4635.0	10676.0	6613.0	1.00
R²_mu[279]	0.64	0.16	0.34	0.96	0.00	0.00	10634.0	10634.0	10524.0	6410.0	1.00

748 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,dataTrace_Rsquared,trace_Rsquared,'R²')

/home/bob/.local/lib/python3.8/site-packages/arviz/data/io_pymc3.py:85: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(
/home/bob/.local/lib/python3.8/site-packages/arviz/data/io_pymc3.py:85: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(
/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(
/home/bob/.local/lib/python3.8/site-packages/arviz/data/io_pymc3.py:85: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(
/home/bob/.local/lib/python3.8/site-packages/arviz/data/io_pymc3.py:85: FutureWarning: Using `from_pymc3` without the model will be deprecated in a future release. Not using the model will return less accurate and less useful results. Make sure you use the model argument or call from_pymc3 within a model context.
  warnings.warn(

with RsquaredModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Rsquared,'R²')

with RsquaredModel as model:
    plotting_lib.pm.energyplot(trace_Rsquared)

Posterior predictive distribution

with RsquaredModel as model:
    posterior_pred_Rsquared = pm.sample_posterior_predictive(trace_Rsquared,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Rsquared,posterior_pred_Rsquared,dataZ["R²_z"].values,'R²')

Compare prior and posterior for model parameters

with RsquaredModel as model:
    pm_data_Rsquared = az.from_pymc3(trace=trace_Rsquared,prior=prior_pred_Rsquared,posterior_predictive=posterior_pred_Rsquared)

arviz.data.io_pymc3 - WARNING - posterior predictive variable R²_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'R²')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Rsquared,'R²')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Rsquared,'R²')

Asfc

with pm.Model() as model:
    AsfcModel = ThreeFactorModel('Asfc',x1,x2,x3,dataZ["Asfc_z"].values)

Verify model settings

AsfcModel.printParams(x1,x2,x3,dataZ["Asfc_z"].values)

The number of levels of the x variables are (2, 11, 140)
The standard deviations used for the beta priors are (1.0422213301584404, 0.6910746316178924, 0.4801483731544936)

pm.model_to_graphviz(AsfcModel)

Check prior choice

with AsfcModel as model:
    prior_pred_Asfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')

Prior choice is as intended: Broad over the data range.

Sampling

with AsfcModel as model:
    trace_Asfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
    #fit_Asfc = pm.fit(random_seed=random_seed)
    #trace_Asfc = fit_Asfc.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [Asfc_b3_dist, Asfc_b2_dist, Asfc_b1_dist, Asfc_b0_dist, Asfc_sigmaY, Asfc_nuY, Asfc_mu_b3, Asfc_mu_b2, Asfc_mu_b1, Asfc_mu_b0, Asfc_sigma3, Asfc_sigma2, Asfc_sigma1, Asfc_sigma0]

100.00% [20000/20000 1:28:26<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 1_000 draw iterations (10_000 + 10_000 draws total) took 5307 seconds.
The number of effective samples is smaller than 10% for some parameters.

with AsfcModel as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('Asfc'), 'wb') as buff:
            pickle.dump({'model': AsfcModel, 'trace': trace_Asfc}, buff)

Check sampling

with AsfcModel as model:
    dataTrace_Asfc = az.from_pymc3(trace=trace_Asfc)

pm.summary(dataTrace_Asfc,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
Asfc_mu_b0	0.06	0.82	-1.51	1.69	0.01	0.01	7447.0	5776.0	7446.0	7394.0	1.0
Asfc_mu_b1[0]	0.03	0.70	-1.33	1.43	0.01	0.01	7109.0	5649.0	7115.0	7132.0	1.0
Asfc_mu_b1[1]	0.01	0.69	-1.32	1.39	0.01	0.01	7552.0	5825.0	7558.0	7622.0	1.0
Asfc_mu_b2[0]	1.02	0.54	-0.09	2.03	0.01	0.01	3499.0	3361.0	3432.0	5699.0	1.0
Asfc_mu_b2[1]	1.17	0.56	0.04	2.24	0.01	0.01	2981.0	2907.0	2952.0	4442.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
Asfc_mu[275]	-1.02	0.03	-1.07	-0.97	0.00	0.00	9456.0	9456.0	9569.0	6577.0	1.0
Asfc_mu[276]	-0.94	0.03	-1.00	-0.88	0.00	0.00	11371.0	11055.0	9957.0	6840.0	1.0
Asfc_mu[277]	-0.97	0.03	-1.02	-0.91	0.00	0.00	11345.0	11109.0	10011.0	6933.0	1.0
Asfc_mu[278]	-1.01	0.02	-1.05	-0.96	0.00	0.00	10217.0	10198.0	10062.0	7580.0	1.0
Asfc_mu[279]	-1.04	0.02	-1.08	-0.99	0.00	0.00	10185.0	10165.0	10023.0	7435.0	1.0

748 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,dataTrace_Asfc,'Asfc')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with AsfcModel as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_Asfc,'Asfc')

with AsfcModel as model:
    plotting_lib.pm.energyplot(trace_Asfc)

Posterior predictive distribution

with AsfcModel as model:
    posterior_pred_Asfc = pm.sample_posterior_predictive(trace_Asfc,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Asfc,posterior_pred_Asfc,dataZ["Asfc_z"].values,'Asfc')

Compare prior and posterior for model parameters

with AsfcModel as model:
    pm_data_Asfc = az.from_pymc3(trace=trace_Asfc,prior=prior_pred_Asfc,posterior_predictive=posterior_pred_Asfc)

arviz.data.io_pymc3 - WARNING - posterior predictive variable Asfc_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_Asfc,'Asfc')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Asfc,'Asfc')

Smfc

with pm.Model() as model:
    SmfcModel = ThreeFactorModel('Smfc',x1,x2,x3,dataZ.Smfc_z.values)

Verify model settings

SmfcModel.printParams(x1,x2,x3,dataZ.Smfc_z.values)

The number of levels of the x variables are (2, 11, 140)
The standard deviations used for the beta priors are (1.3943421938332237, 3.387697874270125, 9.05131085356785)

pm.model_to_graphviz(SmfcModel)

Check prior choice

with SmfcModel as model:
    prior_pred_Smfc = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_Smfc,dataZ.Smfc_z.values,'Smfc')

Prior choice is as intended: Broad over the data range.

Sampling

with SmfcModel as model:
    trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.8,random_seed=random_seed)
    #fit_Smfc = pm.fit(random_seed=random_seed)
    #trace_Smfc = fit_Smfc.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [Smfc_b3_dist, Smfc_b2_dist, Smfc_b1_dist, Smfc_b0_dist, Smfc_sigmaY, Smfc_nuY, Smfc_mu_b3, Smfc_mu_b2, Smfc_mu_b1, Smfc_mu_b0, Smfc_sigma3, Smfc_sigma2, Smfc_sigma1, Smfc_sigma0]

10.79% [2157/20000 17:51<2:27:45 Sampling 10 chains, 0 divergences]

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
   1485             with sampler:
-> 1486                 for draw in sampler:
   1487                     trace = traces[draw.chain - chain]

~/.local/lib/python3.8/site-packages/pymc3/parallel_sampling.py in __iter__(self)
    491         while self._active:
--> 492             draw = ProcessAdapter.recv_draw(self._active)
    493             proc, is_last, draw, tuning, stats, warns = draw

~/.local/lib/python3.8/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout)
    351         pipes = [proc._msg_pipe for proc in processes]
--> 352         ready = multiprocessing.connection.wait(pipes)
    353         if not ready:

/usr/lib/python3.8/multiprocessing/connection.py in wait(object_list, timeout)
    930             while True:
--> 931                 ready = selector.select(timeout)
    932                 if ready:

/usr/lib/python3.8/selectors.py in select(self, timeout)
    414         try:
--> 415             fd_event_list = self._selector.poll(timeout)
    416         except InterruptedError:

KeyboardInterrupt: 

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-85-755e53d074df> in <module>
      1 with SmfcModel as model:
----> 2     trace_Smfc = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=0.8,random_seed=random_seed)
      3     #fit_Smfc = pm.fit(random_seed=random_seed)
      4     #trace_Smfc = fit_Smfc.sample(draws=numSamples)

~/.local/lib/python3.8/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)
    543         _print_step_hierarchy(step)
    544         try:
--> 545             trace = _mp_sample(**sample_args, **parallel_args)
    546         except pickle.PickleError:
    547             _log.warning("Could not pickle model, sampling singlethreaded.")

~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
   1510     except KeyboardInterrupt:
   1511         if discard_tuned_samples:
-> 1512             traces, length = _choose_chains(traces, tune)
   1513         else:
   1514             traces, length = _choose_chains(traces, 0)

~/.local/lib/python3.8/site-packages/pymc3/sampling.py in _choose_chains(traces, tune)
   1528     lengths = [max(0, len(trace) - tune) for trace in traces]
   1529     if not sum(lengths):
-> 1530         raise ValueError("Not enough samples to build a trace.")
   1531 
   1532     idxs = np.argsort(lengths)[::-1]

ValueError: Not enough samples to build a trace.

Analysis stopped here because sampling did not converge. As the plot shows, some data points are very far away from the others, which would require the analysis to be based on more heavy-tailed distributions.

HAsfc9

with pm.Model() as model:
    HAsfc9Model = ThreeFactorModel('HAsfc9',x1,x2,x3,dataZ["HAsfc9_z"].values)

Verify model settings

HAsfc9Model.printParams(x1,x2,x3,dataZ["HAsfc9_z"].values)

The number of levels of the x variables are (2, 11, 140)
The standard deviations used for the beta priors are (1.0299664523718457, 2.0653529113379627, 2.5833026311707217)

pm.model_to_graphviz(HAsfc9Model)

Check prior choice

with HAsfc9Model as model:
    prior_pred_HAsfc9 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')

Prior choice is as intended: Broad over the data range.

Sampling

with HAsfc9Model as model:
    trace_HAsfc9 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
    #fit_HAsfc9 = pm.fit(random_seed=random_seed)
    #trace_HAsfc9 = fit_HAsfc9.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [HAsfc9_b3_dist, HAsfc9_b2_dist, HAsfc9_b1_dist, HAsfc9_b0_dist, HAsfc9_sigmaY, HAsfc9_nuY, HAsfc9_mu_b3, HAsfc9_mu_b2, HAsfc9_mu_b1, HAsfc9_mu_b0, HAsfc9_sigma3, HAsfc9_sigma2, HAsfc9_sigma1, HAsfc9_sigma0]

100.00% [20000/20000 53:56<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 1_000 draw iterations (10_000 + 10_000 draws total) took 3237 seconds.
The rhat statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.
The estimated number of effective samples is smaller than 200 for some parameters.

with HAsfc9Model as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('HAsfc9'), 'wb') as buff:
            pickle.dump({'model': HAsfc9Model, 'trace': trace_HAsfc9}, buff)

Check sampling

with HAsfc9Model as model:
    dataTrace_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9)

pm.summary(dataTrace_HAsfc9,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
HAsfc9_mu_b0	0.02	0.82	-1.57	1.63	0.01	0.01	7718.0	4985.0	7730.0	7316.0	1.00
HAsfc9_mu_b1[0]	0.04	0.69	-1.36	1.36	0.01	0.01	7548.0	5733.0	7518.0	7411.0	1.00
HAsfc9_mu_b1[1]	-0.03	0.70	-1.42	1.30	0.01	0.01	7414.0	5426.0	7407.0	7311.0	1.00
HAsfc9_mu_b2[0]	-0.48	0.52	-1.50	0.56	0.01	0.00	4597.0	4597.0	4582.0	6172.0	1.00
HAsfc9_mu_b2[1]	-0.26	0.52	-1.28	0.73	0.01	0.01	3897.0	3897.0	3880.0	5847.0	1.00
...	...	...	...	...	...	...	...	...	...	...	...
HAsfc9_mu[275]	2.39	0.04	2.31	2.48	0.00	0.00	10050.0	10050.0	10017.0	7890.0	1.00
HAsfc9_mu[276]	2.17	0.36	1.84	3.08	0.01	0.01	761.0	711.0	1589.0	2185.0	1.01
HAsfc9_mu[277]	2.08	0.36	1.75	2.99	0.01	0.01	761.0	709.0	1588.0	2159.0	1.01
HAsfc9_mu[278]	0.27	0.04	0.18	0.36	0.00	0.00	9874.0	9874.0	9900.0	8275.0	1.00
HAsfc9_mu[279]	0.18	0.04	0.10	0.27	0.00	0.00	9767.0	9666.0	9773.0	8522.0	1.00

748 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,dataTrace_HAsfc9,'HAsfc9')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with HAsfc9Model as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc9,'HAsfc9')

with HAsfc9Model as model:
    plotting_lib.pm.energyplot(trace_HAsfc9)

Posterior predictive distribution

with HAsfc9Model as model:
    posterior_pred_HAsfc9 = pm.sample_posterior_predictive(trace_HAsfc9,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc9,posterior_pred_HAsfc9,dataZ["HAsfc9_z"].values,'HAsfc9')

Compare prior and posterior for model parameters

with HAsfc9Model as model:
    pm_data_HAsfc9 = az.from_pymc3(trace=trace_HAsfc9,prior=prior_pred_HAsfc9,posterior_predictive=posterior_pred_HAsfc9)

arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc9_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc9,'HAsfc9')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_HAsfc9,'HAsfc9')

HAsfc81

with pm.Model() as model:
    HAsfc81Model = ThreeFactorModel('HAsfc81',x1,x2,x3,dataZ["HAsfc81_z"].values)

Verify model settings

HAsfc81Model.printParams(x1,x2,x3,dataZ["HAsfc81_z"].values)

The number of levels of the x variables are (2, 11, 140)
The standard deviations used for the beta priors are (1.0035538697180921, 1.674345831114272, 1.8608701769937184)

pm.model_to_graphviz(HAsfc81Model)

Check prior choice

with HAsfc81Model as model:
    prior_pred_HAsfc81 = pm.sample_prior_predictive(samples=numPredSamples,random_seed=random_seed)

plotting_lib.plotPriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')

Prior choice is as intended: Broad over the data range.

Sampling

with HAsfc81Model as model:
    trace_HAsfc81 = pm.sample(numSamples,cores=numCores,tune=numTune,max_treedepth=20, init='auto',target_accept=target_accept,random_seed=random_seed)
    #fit_HAsfc81 = pm.fit(random_seed=random_seed)
    #trace_HAsfc81 = fit_HAsfc81.sample(draws=numSamples)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (10 chains in 10 jobs)
NUTS: [HAsfc81_b3_dist, HAsfc81_b2_dist, HAsfc81_b1_dist, HAsfc81_b0_dist, HAsfc81_sigmaY, HAsfc81_nuY, HAsfc81_mu_b3, HAsfc81_mu_b2, HAsfc81_mu_b1, HAsfc81_mu_b0, HAsfc81_sigma3, HAsfc81_sigma2, HAsfc81_sigma1, HAsfc81_sigma0]

100.00% [20000/20000 1:11:03<00:00 Sampling 10 chains, 0 divergences]

Sampling 10 chains for 1_000 tune and 1_000 draw iterations (10_000 + 10_000 draws total) took 4265 seconds.
The number of effective samples is smaller than 10% for some parameters.

with HAsfc81Model as model:
    if writeOut:
        with open(outPathData + 'model_{}.pkl'.format('HAsfc81'), 'wb') as buff:
            pickle.dump({'model': HAsfc81Model, 'trace': trace_HAsfc81}, buff)

Check sampling

with HAsfc81Model as model:
    dataTrace_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81)

pm.summary(dataTrace_HAsfc81,hdi_prob=0.95).round(2)

	mean	sd	hdi_2.5%	hdi_97.5%	mcse_mean	mcse_sd	ess_mean	ess_sd	ess_bulk	ess_tail	r_hat
HAsfc81_mu_b0	0.03	0.83	-1.60	1.62	0.01	0.01	6328.0	5593.0	6339.0	7044.0	1.0
HAsfc81_mu_b1[0]	0.03	0.70	-1.35	1.38	0.01	0.01	7201.0	5956.0	7198.0	7184.0	1.0
HAsfc81_mu_b1[1]	-0.03	0.69	-1.36	1.35	0.01	0.01	5990.0	5990.0	6007.0	6910.0	1.0
HAsfc81_mu_b2[0]	-0.55	0.54	-1.69	0.47	0.01	0.01	3153.0	3153.0	3150.0	4960.0	1.0
HAsfc81_mu_b2[1]	-0.42	0.54	-1.54	0.58	0.01	0.01	3787.0	3787.0	3774.0	5995.0	1.0
...	...	...	...	...	...	...	...	...	...	...	...
HAsfc81_mu[275]	1.39	0.05	1.29	1.49	0.00	0.00	9779.0	9779.0	9781.0	7388.0	1.0
HAsfc81_mu[276]	2.70	0.26	2.43	3.26	0.01	0.00	1242.0	1169.0	2138.0	4486.0	1.0
HAsfc81_mu[277]	2.62	0.26	2.36	3.18	0.01	0.00	1241.0	1166.0	2111.0	4217.0	1.0
HAsfc81_mu[278]	0.65	0.10	0.48	0.82	0.00	0.00	7362.0	7321.0	7854.0	8632.0	1.0
HAsfc81_mu[279]	0.57	0.10	0.40	0.74	0.00	0.00	7379.0	7313.0	7925.0	8249.0	1.0

748 rows × 11 columns

plotting_lib.plotDiagnostics(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,dataTrace_HAsfc81,'HAsfc81')

/home/bob/.local/lib/python3.8/site-packages/arviz/plots/backends/matplotlib/pairplot.py:212: UserWarning: rcParams['plot.max_subplots'] (40) is smaller than the number of resulting pair plots with these variables, generating only a 8x8 grid
  warnings.warn(

with HAsfc81Model as model:
    plotting_lib.plotTracesB(widthInch,heigthInch,dpi,writeOut,outPathPlots,trace_HAsfc81,'HAsfc81')

with HAsfc81Model as model:
    plotting_lib.pm.energyplot(trace_HAsfc81)

Posterior predictive distribution

with HAsfc81Model as model:
    posterior_pred_HAsfc81 = pm.sample_posterior_predictive(trace_HAsfc81,samples=numPredSamples,random_seed=random_seed)

/home/bob/.local/lib/python3.8/site-packages/pymc3/sampling.py:1707: UserWarning: samples parameter is smaller than nchains times ndraws, some draws and/or chains may not be represented in the returned posterior predictive sample
  warnings.warn(

100.00% [2000/2000 00:02<00:00]

plotting_lib.plotPriorPosteriorPredictive(widthInch,heigthInch,dpi,writeOut,outPathPlots,df,dictMeanStd,prior_pred_HAsfc81,posterior_pred_HAsfc81,dataZ["HAsfc81_z"].values,'HAsfc81')

Compare prior and posterior for model parameters

with HAsfc81Model as model:
    pm_data_HAsfc81 = az.from_pymc3(trace=trace_HAsfc81,prior=prior_pred_HAsfc81,posterior_predictive=posterior_pred_HAsfc81)

arviz.data.io_pymc3 - WARNING - posterior predictive variable HAsfc81_y's shape not compatible with number of chains and draws. This can mean that some draws or even whole chains are not represented.

plotting_lib.plotPriorPosteriorB(widthInch,heigthInch,dpi,sizes,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')

Posterior and contrasts

plotting_lib.plotPosterior(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,pm_data_HAsfc81,'HAsfc81')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_HAsfc81,'HAsfc81')

Summary

The contrast plots between the two software packages are shown below again for each variable except Smfc

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_epLsar,'epLsar')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Rsquared,'R²')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_Asfc,'Asfc')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_HAsfc9,'HAsfc9')

plotting_lib.plotContrast(widthInch,heigthInch,dpi,writeOut,outPathPlots,dictMeanStd,x1contrast_dict,trace_HAsfc81,'HAsfc81')

Write out

!jupyter nbconvert --to html Statistical_Model_ThreeFactor.ipynb

[NbConvertApp] Converting notebook Statistical_Model_ThreeFactor.ipynb to html
[NbConvertApp] Writing 9742523 bytes to Statistical_Model_ThreeFactor.html

!jupyter nbconvert --to markdown Statistical_Model_ThreeFactor.ipynb

[NbConvertApp] Converting notebook Statistical_Model_ThreeFactor.ipynb to markdown
[NbConvertApp] Support files will be in Statistical_Model_ThreeFactor_files/
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Making directory Statistical_Model_ThreeFactor_files
[NbConvertApp] Writing 81991 bytes to Statistical_Model_ThreeFactor.md